<html>
<head>
<link rel="stylesheet" href="../blueprint/screen.css" type="text/css" media="screen, projection">

<title>Raster Manipulation In Python</title>

<link rel="stylesheet" href="../css/styles.css" type="text/css" media="screen,print">

</head>
<body>
<div class="container">
<h1>Raster Manipulation in Python</h1>

<p>
The <a href="index.html">rasterlang</a> plugin provides a language for 
raster manipulation, but if you want to do something more complex then
it provides a couple of useful functions for working with rasters in 
Python.
</p>

<h2>Read a Raster into an array</h2>
<p>
Start up the Python console and load some rasters. Then in the console do:
<pre>
layer1 = iface.mapCanvas().layer(0)
layer2 = iface.mapCanvas().layer(1)
layer1.name()
layer2.name()
</pre>
This should show the names of your raster layers. Next we convert them to arrays:
<pre>
from rasterlang.layers import layerAsArray
a1 = layerAsArray(layer1)
a2 = layerAsArray(layer2)
a1.shape
a2.shape
</pre>
The shape property tells you the number of bands, rows, and columns of your raster. If the raster has only one band then the shape only has two values.
</p>

<h2>Computing with arrays</h2>
<p>
For information on calculating with arrays, see the <a aref="http://numpy.scipy.org">Numpy Documentation</a>. You can do <code>import numpy</code> in the console and use the numpy functions. Numpy also makes life easy by defining the usual arithmetic operators on arrays. For example, we can add our rasters:
<pre>
rsum = a1 + a2
</pre>
Note that numpy will raise an exception if the arrays aren't
conformable, which basically means they have matching numbers of rows and
columns. It's all explained in the numpy docs. You can add a single band raster 
to a multiple band raster if they have the same number of rows and columns -
the single band is added to all of the multiple bands.
</p>

<h2>Saving Arrays as GeoTiffs</h2>
<p>
You can save an array to a GeoTiff using another function from rasterlang's array
module. Note that arrays do not have any geographical location, so you need to 
supply an 'extent' - this is a list of length 4 giving <code>[xmin,ymin,xmax,ymax]</code>. It is most easily made from an existing layer using the <code>.extent()</code> method. 

<pre>
from rasterlang.layers import writeGeoTiff
e = layer1.extent()
extent = [e.xMinimum(),e.yMinimum(),e.xMaximum(),e.yMaximum()]
writeGeoTiff(rsum, extent, "rsum.tiff")
</pre>

You can then load this layer into Qgis using the menu options or a bit more Python code:

<pre>
newLayer = QgsRasterLayer("rsum.tiff","rsum")
QgsMapLayerRegistry.instance().addMapLayer(newLayer)
</pre>

</p>

<h2>A Note on Numpy and Numeric</h2>
<p>

Originally there was the <code>numeric</code> package for Python array handling. Then this developed into <code>numpy</code>, which is closely compatible but not exactly interchangeable with <code>numeric</code>. The gdal package for raster data file handling is transitioning from <code>numeric</code> to <code>numpy</code>, so if you have an older gdal then you may have some problems...
</p>
<p>
What happens is that you read the layer as an array and get back a <code>numeric</code> array from gdal. Then you manipulate it with <code>numpy</code> functions and it becomes a <code>numpy</code> array (<code>numpy</code> has done the conversion). Then when you go to save the GeoTiff, you are passing a <code>numpy</code> array to a gdal that expects a <code>numeric</code> array, and it fails.
</p>
<p>
The solution is to convert the resulting array back to the right type. If you need to turn a <code>numpy</code> array into a <code>numeric</code> one, then do:
<pre>
import Numeric
aNumer = Numeric.array(anumpy)
</pre>
The rasterlang plugin works this out for itself, but if you are doing this in 
Python you'll have to do the conversion if necessary. Once every gdal is using purely <code>numpy</code> then this won't be needed.
</p>

<h2>Credits</h2>
<p>
Written by Barry Rowlingson &lt;b.rowlingson@lancaster.ac.uk&gt;
</p>
</div>
</body>
</html>
